! Project TWS Matrixexponentiële
!
!
! MODULE METHODS
! 
! Deze module bevat de verschillende algoritmen
!
! @author Roel Matthysen
! @version 2.0
module methods
    
    use fftw
    use IO
    implicit none
    private

    integer, parameter :: doubleprecision = kind(0.d0)
    real(doubleprecision), parameter :: tol = 10.0d0**(-10)
    public :: taylorexp,eig,eigqr,taylorexpc,eigc, eigqrc
    
    contains
            
        ! TAYLOREXP
        !
        ! Deze functie berekent de exponentiële van een volle matrix A met de taylor methode.
        ! 
        ! @param    A                    De volle invoermatrix A
        !
        ! @param    verbose(optional)    Dit argument bepaalt of debugging-informatie wordt afgedrukt.
        !
        ! @result   RES                  Een tweedimensionale doubleprecision matrix met de exponentiële van A.
        !
        ! @since    1.0

        function taylorexp(A,verbose) result(RES)
            real(doubleprecision), dimension(:,:), intent(in) :: A
            logical, optional ,intent(in) :: verbose
            real(doubleprecision) :: RES(size(A,1),size(A,2)),add(size(A,1),size(A,2))
            integer :: i
            ! Druk Matrix A af
            if(verbose) print *, size(A,1)
            if(verbose) print *, size(A,2)
            if(verbose) call printmatrix(A,'A')
            ! Add en RES worden op nul geïnitialiseerd.
            add = 0.0
            RES = 0.0
            ! Add is een eenheidsmatrix
            do i = 1, size(add,1)
                add(i,i) = 1
            enddo
            ! RES wordt berekend zoals in de referentie-implementatie
            do i = 1,1000
                    RES = RES+add;
                    add = matmul(add,A)/i;
                    if(maxval(abs(add/RES))<=tol) then
                        if(verbose) call printmatrix(add,'add')
                        exit
                    endif
            enddo
        end function
        
        
        ! EIG
        !
        ! Deze functie berekent de exponentiële van een volle matrix A met de taylor eigenwaarden methode.
        ! 
        ! @param     A                    De volle invoermatrix A
        !
        ! @param     verbose(optional)    Dit argument bepaalt of debugging-informatie wordt afgedrukt.
        !
        ! @result    RES                  Een tweedimensionale doubleprecision matrix met de exponentiële van A.
        !
        ! @since     1.0
        function eig(A,verbose) result(RES)
            real(doubleprecision), dimension(:,:), intent(in) :: A
            logical, optional, intent(in) :: verbose
            integer :: j,i, ipvt(size(A,1)),N, worksize
            real(doubleprecision) :: RES(size(A,1),size(A,1))
            complex(doubleprecision) :: w(size(A,1)),vl(size(A,1),size(A,1))
            complex(doubleprecision) :: vr(size(A,1),size(A,1))
            N = size(A,1)
            vl=0.0
            ! Eigenwaarden berekenen
            call zgeeveigen(A,w,vr, worksize)
            if(verbose) call printmatrix(vr,'eigenvectoren')
            if(verbose) call printmatrix(w,'eigs')
            w = exp(w)
            if(verbose) call printmatrix(w,'eigsexp')
            ! Exponentiëlen op de diagonaal
            do i=1,N
                vl(i,i)=w(i)    
            enddo
            if(verbose) call printmatrix(vl,'expdiag')
            vl = matmul(vr,vl)
            vl = transpose(vl)
            if(verbose) call printmatrix(vl,'vl')
            vr = transpose(vr)
            ! Stelsel oplossen
            call ZGESV(N,N,vr,N,ipvt,vl,N,j)
            if(.NOT. j==0) print *, 'Error in ZGESV'
            RES = dble(transpose(vl))
            if(verbose) call printmatrix(RES,'Resultaat')
        end function
        
        ! EIGQR
        !
        ! Deze functie berekent de exponentiële van een volle matrix A met de taylor eigenwaarden methode + QR-o
        ! ontbinding.
        ! 
        ! @param     A                    De volle invoermatrix A
        !
        ! @param     verbose(optional)    Dit argument bepaalt of debugging-informatie wordt afgedrukt.
        !
        ! @result    RES                  Een tweedimensionale doubleprecision matrix met de exponentiële van A.
        !
        ! @since     1.0
        function eigqr(A,verbose) result(RES)
            real(doubleprecision), dimension(:,:), intent(in) :: A
            logical, optional, intent(in) :: verbose
            integer :: j,i,N, worksize
            real(doubleprecision) :: RES(size(A,1),size(A,1))
            complex(doubleprecision) :: tau(size(A,1))
            complex(doubleprecision), allocatable :: workqr(:)
            complex(doubleprecision) :: w(size(A,1)),vl(size(A,1),size(A,1))
            complex(doubleprecision) :: vr(size(A,1),size(A,1))
            complex(doubleprecision) :: eigv(size(A,1),size(A,1))
            N = size(A,1)
            vl = 0.0
            eigv = 0.0
            worksize=0
            w = 0.0
            ! Eigenwaarden berekenen
            call zgeeveigen(A,w,eigv,worksize)
            vr = eigv
            allocate(workqr(worksize))
            ! QR-ontbinding
            call ZGEQRF(N,N,vr,N,tau,workqr,worksize,j)
            if(.NOT. j==0) print *, 'Error in ZGEQRF'
            if(verbose) then
                call printmatrix(vr,'Result of QR')
                call printmatrix(w,'eigs')
                call printmatrix(eigv,'eigv')
            endif
            w = exp(w)
            ! Exponentiëlen op de diagonaal
            do i=1,N
                vl(i,i)=w(i)    
            enddo
            ! vl blijft hetzelfde
            vl = matmul(eigv,vl)
            if(verbose) call printmatrix(vl,'V*exp(eigs)')
            eigv = vr
            ! Q expliciet opstellen
            call ZUNGQR(N,N,N,eigv,N,tau,workqr,worksize,j)
            if(.NOT. j==0) print *, 'Error in ZUNGQR'
            if(verbose) call printmatrix(eigv,'Q')
            eigv = transpose(conjg(eigv))
            ! RY=Q' triangulair systeem oplossen
            call ZTRTRS('U','N','N',N,N,vr,N,eigv,N,j)
            if(.NOT. j==0) print *, 'Error in ZTRTRS'
            RES = real(matmul(vl,eigv))
            if(verbose) then
                call printmatrix(eigv,'Y')
                call printmatrix(matmul(vl,eigv),'imagresult')
                call printmatrix(RES,'final result')
            endif
            deallocate(workqr)
        end function
        
        
        ! TAYLOREXPC
        !
        ! Deze functie berekent de exponentiële van een circulante matrix met de taylor methode.
        ! 
        ! @param     A                    De eerste kolom van de circulante invoermatrix
        !
        ! @param     verbose(optional)    Dit argument bepaalt of debugging-informatie wordt afgedrukt.
        !
        ! @result    RES                  Een ééndimensionale doubleprecision matrix met de eerste kolom van  de 
        !                                 exponentiële van de circulante invoermatrix.
        !
        ! @since     2.0
        function taylorexpc(Ac,verbose) result(RES)
            real(doubleprecision), dimension(:), intent(in) :: Ac
            logical, optional, intent(in) :: verbose
            real(doubleprecision) :: RES(size(Ac,1))
            complex(doubleprecision) :: SUM(size(Ac,1)),add(size(Ac,1)), c(size(Ac,1))
            integer :: i
            ! Druk Matrix A af
            if(verbose) call printmatrix(Ac,'Ac')
            ! Add en RES worden op nul geïnitialiseerd.
            add = 0.0
            SUM = 0.0
            ! Add is een eenheidsmatrix
            do i = 1, size(add,1)
                add(i) = 1
            enddo
            c=fft(Ac)
            if(verbose) call printmatrix(c,'c')
            !RES wordt berekend zoals in de referentie-implementatie
            do i = 1,1000
                    SUM = SUM+add;
                    add = add*c/i;
                    if(maxval(abs(add/SUM))<=tol) then
                        if(verbose) call printmatrix(add,'add')
                        exit
                    endif
            enddo
            RES=ifft(SUM)
        end function
        
        ! EIGC
        !
        ! Deze functie berekent de exponentiële van een circulante matrix met de eigenwaardenmethode.
        ! 
        ! @param     A                    De eerste kolom van de circulante invoermatrix
        !
        ! @param     verbose(optional)    Dit argument bepaalt of debugging-informatie wordt afgedrukt.
        !
        ! @result    RES                  Een ééndimensionale doubleprecision matrix met de eerste kolom van de      
        !                                 exponentiële van de circulante matrix.
        !
        ! @since     2.0
        function eigc(Ac,verbose) result(RES)
            real(doubleprecision), dimension(:), intent(in) :: Ac
            logical, optional, intent(in) :: verbose
            real(doubleprecision) :: RES(size(Ac,1))
            complex(doubleprecision) :: c(size(Ac,1))
            c=fft(Ac)
            if(verbose) call printmatrix(c,'FFT Transformatie van inputmatrix A')
            c=exp(c)
            RES=ifft(c)
        end function
        
        ! EIGQRC
        !
        ! Deze functie berekent de exponentiële van een circulante matrix met de eigenwaardenmethode + QR-ontbinding.
        ! 
        ! @param     A                    De eerste kolom van de circulante invoermatrix
        !
        ! @param     verbose(optional)    Dit argument bepaalt of debugging-informatie wordt afgedrukt.
        !
        ! @result    RES                  Een ééndimensionale doubleprecision matrix met de eerste kolom van de      
        !                                 exponentiële van de circulante matrix.
        !
        ! @since     2.0
        function eigqrc(Ac,verbose) result(RES)
            real(doubleprecision), dimension(:), intent(in) :: Ac
            logical, optional,intent(in) :: verbose
            real(doubleprecision) :: RES(size(Ac,1))
            complex(doubleprecision) :: c(size(Ac,1))
            c=fft(Ac)
            if(verbose) call printmatrix(c,'FFT Transformatie van inputmatrix A')
            c=exp(c)
            RES=ifft(c)
        end function
    
        ! ZGEEVEIGEN
        !
        ! Deze functie berekent de eigenwaarden en eigenvectoren van een matrix A
        ! 
        ! @param     A                    De volle matrix
        !
        ! @param     w                    bevat na afloop de eigenwaarden
        !
        ! @param     vr                   bevat na afloop de eigenvectoren
        !
        ! @param     worksize             bevat na afloop de optimale worksize
        !
        ! @since     2.0        
        subroutine zgeeveigen(A,w,vr, worksize)
            real(doubleprecision), dimension(:,:) :: A
            real(doubleprecision), allocatable :: Acopy(:,:)
            complex(doubleprecision) :: w(:)
            complex(doubleprecision) :: vr(:,:)
            real(doubleprecision) :: dvr(size(A,1),size(A,1)), dvl(size(A,1),size(A,1)), wr(size(A,1))
            real(doubleprecision) :: wi(size(A,1))
            real(doubleprecision), allocatable :: dwork(:)
            integer :: i,j, N, worksize
            logical :: eigvtodo
            eigvtodo = .true.
            N = size(A,1)
            ! Deze kopie is nodig om A onveranderd te laten. Aan DGEEV kan (A) meegegeven worden,
            ! maar ifort geeft hier warnings voor.
            allocate(Acopy(N,N))
            Acopy = A
            allocate(dwork(1))
            call DGEEV('V','V',N,Acopy,N,wr,wi,dvl,N,dvr,N,dwork,-1,j)
            worksize=aint(real(dwork(1)));
            deallocate(dwork)
            allocate(dwork(worksize))
            if(.NOT. j==0) print *, 'Error in DGEEV'
            call DGEEV('V','V',N,Acopy,N,wr,wi,dvl,N,dvr,N,dwork,worksize,j)
            do i=1,N
                if(wi(i)==0) then
                    w(i) = cmplx(wr(i),0.d0,kind=doubleprecision)
                    do j=1,N
                        vr(j,i) = cmplx(dvr(j,i),0.d0,kind=doubleprecision)
                    enddo
                else
                    w(i) = cmplx(wr(i),wi(i),kind=doubleprecision)
                    if(eigvtodo) then
                        do j=1,N
                            vr(j,i) = cmplx(dvr(j,i),dvr(j,i+1),kind=doubleprecision)
                            vr(j,i+1) = cmplx(dvr(j,i),-1.d0*dvr(j,i+1),kind=doubleprecision)
                        enddo
                        eigvtodo = .false.
                    else
                        eigvtodo = .true.
                    endif
                endif
            enddo
            deallocate(dwork)
            deallocate(Acopy)
        end subroutine
        
        
        
end module
